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The density functional theory (DFT) interaction energy of a dimer is rigorously derived from the 
monomer densities. To this end, the supermolecular energy bifunctional is formulated in terms of 
mutually orthogonal sets of orbitals of the constituent monomers. The orthogonality condition is 
preserved in the solution of the Kohn-Sham equations through the Pauli blockade method. Numer- 
ical implementation of the method provides interaction energies which agree with those obtained 
from standard supermolecular calculations within less than 0.1 % error for three example function- 
als: Slater-Dirac, PBEO and B3LYP, and for two model van der Waals dimers: Ne2 and (C2H4)2, 
and two model H-bond complexes: (HF)2 and (NH3)2. 



I. INTRODUCTION 

Deriving the dimer interaction energy via mutual po- 
larization of constituent monomers is important both 
from the fundamental perspective and from a practical 
point of view. In particular, it may aid the understand- 
ing how the non-covalent systems are described in the 
density functional theory which is one of the most prob- 
lematic issues of the electronic structure theory. The 
major problem of DFT as applied to the van der Waals 
systems is a wrong description of dispersion forces. "'^^'^ 
Surprisingly enough, little has been done to better un- 
derstand the performance of supermolecular interaction 
energy in the framework of DFT. For the Hartree-Fock 
(HF) interaction energy, such an approach has been 
pioneered by Morokuma"' in the 1970s, and a decade 
later, inspired by the work of Sadlej^, rigorously derived 
by Gutowski and Piela^ (see also Ref.^). The pertur- 
bation approach within the symmetry-adapted pertur- 
bation theory (SAPT) formalism was also extensively 
exploited in this context. ^'^ In the age of DFT, it is 
highly desirable to develop such an approach also for 
the density functional formalism. Approximate DFT 
treatments have already been advanced by Cortona^^ 
and Wesolowski and Warshel^^; see also recent energy 
decomposition schemes proposed in Refs. Cybulski and 
Seversen^^, Reinhardt et al.^'^, Su and Li^'*, and Refs. 
therein, as well as the density functional formulation of 
SAPT.i'5'16 

The goal of this work is to derive rigorously the super- 
molecular density functional theory (DFT) interaction 
energy via the mutual polarization of the monomer densi- 
ties. To this end the supermolecular (dimer) energy func- 
tional is expressed in terms of mutually orthogonalized 
sets of the Kohn-Sham (KS) orbitals of the constituent 
monomers. The coupled KS equations are next solved 
iteratively, by using the Pauli blockade (PB) technique 
of Gutowski and Piela^. The correctness of the deriva- 
tion is demonstrated by comparing the interaction energy 
calculated from the equations developed here and the su- 



permolecular interaction energies. The DFT approxima- 
tion to the Heitler London interaction energy, based on 
the decomposition of the interaction energy introduced 
in this paper, is also discussed. 



II. THEORY 

In this paper we consider the interaction between two 
closed-shell systems, however, the generalization for high- 
spin open-shell systems and clusters is straightforward. 
The supermolecular interaction energy in terms of DFT 
can be defined as the difference between the total energies 
of the dimer AB and the individual monomers A and B, 
separated to infinity: 



TTiDFT _ pDFT pDFT pDFT 
^int — ^AB ^ ^ 



(1) 



It was demonstrated by Gutowski and Piela^ that the HF 
supermolecular interaction energy may be exactly recov- 
ered by solving the HF equations for monomers in the 
presence of the external perturbation, consisting of the 
electrostatic potential and the non-local exchange poten- 
tial generated by the second monomer. They have also 
proposed a convenient computational scheme for the PB 
method in terms of mutually orthogonalized A and B 
orbitals with the penalty operator forcing the orthogo- 
nality of monomers' occupied orbitals. In this section we 
derive an analogous formalism in terms of appropriately 
modified KS equations and monomer densities. 

We begin with KS equations for the isolated monomers 
which yield the starting orbitals {0"}^^^ and {b^}f.pg- 
The orbitals of monomer A are the solutions of the fol- 
lowing eigen equation: 



/^^■°(r)a?(r) = ei,a°(r). 



(2) 



and the analogous equations hold for the monomer B. 
The KS operator of Eq. (2) is written as 



2 " 



(3) 
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where the monomer A nuclear potential is 



E 



(4) 



with .yVp^ being the number of monomer A nuclei, each 
described by its position and charge Za- Its coulom- 
bic potential reads 



iA(r) 



\v-r\ 

The total energy of monomer A can be written as 
Ea [pI] = [pI] + Vr [p1] + J [pI] + 

rnr 



(5) 



(6) 



The functional (6) comprises the non-interacting kinetic 
energy: 



(7) 



with A being the set of the indices of occupied orbitals 
of the monomer A, nuclear-electron attraction energy: 



coulombic energy 



<°(r)p^(r)d3r, 



r\2 



(8) 



(9) 



exchange-correlation (xc) energy 

[pa] = F-\p\^^): {VrP^(r); . . . }) d\ (10) 

which is evaluated through the numerical integration of 
the F^'^ integrand on a grid of points around monomer A, 
and the nuclear-nuclear repulsion term 



E E 



ZnZ, 



(11) 



which is constant for a fixed geometry. The density of 
monomer A is 



pi(r) = 2^|a°(r) 



(12) 



Similar expressions can be written for monomer B. 

The original, isolated-monomer orbital sets {a^jjg^ 
and are not mutually orthogonal. To proceed, it 

is also important to introduce the set of orthonormalized 
orbitals which are obtained by using Lowdin symmetric 
orthonormalization.^^ The quantities expressed in the or- 
thonormalized orbitals are henceforth marked with tilde. 



One should remember that the orthonormalization leaves 
the total density of the dimer unchanged. However, it 
does change the monomer densities into the densities de- 
formed by the presence of the interacting partner. 

In the PB method the zeroth-order wavefunction of 
the dimer is the wavefunction of the system in the ab- 
sence of molecular interaction. It is constructed from the 
antisymmetrized product of the orthogonalized occupied 
orbitals of the monomers A and B. In case of KS equa- 
tions for dimer the KS determinant can be constructed 
in the same manner as: 



^AB = ^i'li^l, 



(13) 



where V'a and t/'b are KS determinants of monomers A 
and B, respectively. Since the determinants are con- 
structed from orthonormalized orbitals, the V'ab "^ot:- 
malized. 

It can be easily shown that the zeroth-order density of 
the system can be simply written as a sum of monomer 
densities expressed in terms of orthonormalized orbitals. 



Pab^Pab^~Pa + pI- 



(14) 



Note that (14) does not hold for the densities obtained 
from nonorthonormal orbitals of the dimer, i.e. ^ 



pI 



p'b- 



If the interaction between monomers is switched on we 
assume that the KS determinant of the dimer can be writ- 
ten as the antisymmetrized product of two determinants 
for both the monomers: 



AB 



(15) 



and hence the dimer density fulfills the additivity con- 
dition (14). Owing to (14) and using (6), the energy 
functional for the system corresponding to (15) is 

-Bab[pab] = EabIpab] EabIpa + Pb] = 
= T'[pA + Pb] + VXb[pa + Pb] + J[pa + Pb] + 
+ E-^[pa + pb] + VI^. (16) 

Now we rewrite the functional (IG) extracting the 
monomer contributions to the dimer energy through a 
careful inspection of the terms in (16). It is clear from (7) 
that the non-interacting kinetic energy functional is lin- 
ear, 

TnPA + PB]=T'[pA]+T'[pB], (17) 
the nuclear-electron attraction may be separated as 

^a"b[pa + Pb] = 

(<°(r) + v^%r)) (pA(r) + Mr)) d'r = 

= vnpA] + vi''[pB]+ 

+ f v^^ir)pA{r)d'r+ f <(r)pB(r) d^r, (18) 
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and the coulombic term may be decomposed according 
to 



-1 ^3„^ ^3 



[pK^yi) + PB(ri)^ (pA(r2) + PB(r2 

JR3 -/RS ''12 



(19) 



However, the explicit analytical form of the xc func- 
tional is unknown and its approximations depend on the 
functional used. Thus, we introduce the xc energy non- 
additivity, Ai^xc^ 

^E^,\pt^^M^E^'\(>K^PB\-E^'\pP^\-E'"'\pB\. (20) 

It is worthwhile to note that the present formulation nei- 
ther separates nor approximates any of the kinetic non- 
additivity terms appearing in the method of Wesolowski 
and Warshel^"'^. These terms are implicitly and exactly 
included in the term (20) and thus are automatically ac- 
counted for in a consistent manner for any functional. 
Although the expression (16) is a functional of a single 
density, we now make use of (14) and treat the system en- 
ergy as a bifunctional depending on both monomer den- 
sities: 



-Eab[pa + Pb] = £^ab[pa;pb] 



(21) 



Thus, in our search for the ground-state dimer energy, we 
will minimize, with respect to pA and /5b, the bifunctional 
of the form: 

EkS,\PK\PYi\ = 

= t^[pa] -I- vr{PK\ + j{pk\ + E-^{p^\ + vr+ 

+ £^int[PA; Pb] 

= EpXpk\ + £;b[pb] + ^int[PA; pb], (22) 



where 



-EintiPA; Pb] = 

= / «g^(r)/5A(r)d3r+ / vtivYp^ir) d^v^ 



PA(ri)pB(r2) ,3^ , ,,nn 



'"12 



A£;xc[pa;pb] = 

£^cist[pA; pb] + A£'xc[pa; pb] 



(23) 



In the above equation, V-^ is intermonomer nuclear- 
nuclear repulsion energy. However, for the density addi- 
tivity condition (14) to hold, all orbitals must be kept mu- 
tually orthogonal, and the orthogonality also ensures that 
the intersystem Pauli exclusion principle is fulfilled. To 



this end, we perform the variational optimization in two 
steps, using the Pauli blockade (PB) method of Gutowski 
and Piela^: first, the bifunctional extremal search is per- 
formed without the imposition of the intermonomer or- 
thogonality constraint, and secondly, the penalty oper- 
ator is added in the resulting iterative scheme. The 
penalty operator for monomer A reads 



^A = X! I"* 



(24) 



and it is obvious that its action on monomer B's orbitals 
annihilates them once the orbitals are orthogonal. Now 
we turn to the first step: to find a bifunctional mini- 
mum, we calculate the variational derivative of (22) with 
respect to pp^: 

'^-5'ab[pa;Pb] _ 
^PA(r) 



1 



Ar + <(r)+jA(r)-H^r(r)+ 



+ i;J^°(r)+jB(r) + A^r(r) = 
= /F(r) + A<(r)+^^'^*(r), 

where the electrostatic potential is 

^B^*(r)=«B(r)+jB(r) 

and the non-additivity xc operator reads 

A«r(r) = 
^ 6^E^-\pk;M ^ 
SpA{r) 

SE^'^Ipa + pb] SE-<=[pa] 



(25) 
(26) 



'5pa(i-) 

SE^'^Ipab] 



^PA{r) 



(5pAB(r) 

= vTBir)-vTir). (27) 
Hence, the Euler equation for the bifunctional (22) is 

1 



PA 



-Ar + i)f(r) 



(28) 



with 

^■.eff i".clst /„\ I --XC /„\ I -cist 



vf{r) = iT\r) + vTir) + vT{r) + iTiv) + A^r(r), 

(29) 

and Pa being the Lagrange multiplier for the constraint: 



/ PA{r)dh = 



(30) 



The minimization of (22) with respect to pB proceeds in 
an analogous way. Finally, the orbitals minimizing the 
functional (22) are determined to satisfy 



^/KS(r) + A«X^(r) + 5^i''*(r)ja,(r) = eA,.a.(r) 
7f S(r) + A«g-(r) + it\r))h{v) = eji/bk{v) 



(31) 



4 



In the second step of the PB procedure, the iterative 
process of solving Eqs. (31) with the aid of the penalty 
operator is formulated. For monomer A the nth iterative 
step reads 



' ?KS[n-l] 



fa 



-clst[n— 1] 



(32) 



and its equivalent for monomer B is obtained through 
the interchange of the A and B subscripts in (32). 77 > 
is a parameter not affecting the final solutions. The or- 
bitals obtained in (32) are orthogonalized, yielding an 

orthonornial l^lai"'! ;|^{r'| j" set. The interac- 
tion energy at the nth iteration is obtained upon the 
insertion of the densities calculated with these orbitals 
into (22) and subtracting the unperturbed monomer en- 
ergies: 



pPB[n] _ p 



~n ~\n\ 



~[n\ ~[7i 

Pa ;Pb 



Pa ;pb 



(33) 



In the above equation, the A monomer deformation is 



Pa 



Ea [pI] , 



(34) 



and analogously for monomer B. Upon reaching self- 
consistency, the energy (33) is equal to the supermolec- 
ular DFT interaction energy of (1). The computational 
cost of our approach is essentially the same as that of the 
standard KS calculations. 

Since the iterative process (32) starts with the KS or- 
bitals of the isolated monomers, the zero-iteration in- 
teraction energy may be viewed as an analog of the 
well-known HF-based Heitler-London interaction energy. 
Specifically, we define the DFT-based HL interaction en- 
ergy as 



E, 



HL 



E: 



PB[0] 



Eab [pI;p%\^Ea[pI]-Eb[p'b] 



(35) 

This definition is equivalent to that proposed by Cy- 
bulski and Seversen^^. The difference between the self- 
consistent interaction energy and the HL interaction en- 
ergy. 



pPB 



E, 



PB 



E, 



HL 



(36) 



is referred to as the deformation energy. 

It should be stressed here that both and -E™, as 
defined by the above equations, are uniquely defined, and 
are independent of the orthogonalization procedure, and 
may be interpreted in terms of SAPT. The E^ is the HL 
energy arising between the unperturbed, non-orthogonal 
monomers. At the HF level of theory, it includes the in- 
termolecular electrostatic and exchange energies. In the 
DFT case, depending on a particular functional, it may 



also contain some obscure residual inter-monomer elec- 
tron correlation terms that are related to the dispersion 
effect. This is because the interaction operator is, in gen- 
eral, the exchange-correlation operator rather than the 
exact-exchange one only, and the correlation is basically 
of a local type. 

The E^^ term represents the deformation effect with 
respect to the non-orthogonal isolated monomers. First, 
it includes both the induction effects and the CT effects 
that are related to the induction and exchange-induction 
energies as defined by SAPT except that it is obtained 
iteratively through the infinite order rather than pertur- 
batively through the finite order. Second, if an exact 
exchange-correlation operator were used the dispersion 
energy would be included in -E™. As for now, the major- 
ity of existing functionals fail to to account for dispersion 
and at the same time they are not entirely dispersion- 
free. Consequently, for such functionals, E^^ contains 
some residual dispersion terms as well. 



III. NUMERICAL RESULTS 
A. Computational Details 

The method described in Sec. II has been coded 
within the Molpro program suite. Numerical calcu- 
lations were carried out for three different function- 
als: Slater-Dirac^^ (henceforth termed DIRAC), PBEO,^" 
and B3LYP,^^^^^ and for model systems: two van der 
Waals complexes (Ne2 and the ethylene dimer), and 
two hydrogen-bonded dimers [(HF)2 and (NH3)2]. For 
comparison, we also present the results for the stan- 
dard Hartree-Fock method (HF). The distance between 
Ne atoms in Ne2 was set to 6 oq. The geometries for 
(NH3)2 and (C2H4)2 were taken from Ref.'^'^ and from 
Ref.'^* for (HF)2. The numerical procedure depends on 
the following parameters: energy threshold, i.e. the min- 
imum difference in interaction energies from successive 
iterations for which the iterations are continued, the grid 
threshold, i.e. the accuracy with which the Slater-Dirac 
functional can be integrated using the grid as compared 
to its analytical integral, and rj parameter of Eq. (32). 
77 = 10^ was used for all calculations. However, we stress 
once more that rj does not affect the final solutions, only 
the convergence. 

In Table I we present the numerical values of the com- 
ponents of Eq. (33) together with the DFT-based HL 
interaction energy (35), the deformation (36), and the 
relative difference between the interaction energies ob- 
tained in a bifunctional and supermolecular approaches. 



6E„ 



rpPB _ pDFT 

'"^^dft"^* • 100%. 



(37) 



In Fig. 1 we present the dependence of the relative dif- 
ference (37) on the grid threshold. The energy thresh- 
old was set to 10~^ mH and it was kept at that value 
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TABLE I: Interaction energies and their contributions for the bifunctional approach and its comparison with the DFT super- 
molecular energies. All values in mH. The numbers in parentheses denote powers of ten. 



System 


Functional AEa 


AEb 


-Bolst 


ASxc 




tttPB 
^dcf 


TJlPB 

Eint 


SEint/% 




DIRAC 


0.658 


0.658 


-0.953 


-0.625 


-0.262 


-0.0249 


-0.287 


-5.54(-5) 


Nea 


PBEO 


0.416 


0.416 


-0.587 


-0.343 


-0.098 


-0.00862 


-0.107 


-3.74(-6) 


{R = 6 ao) 


B3LYP 


0.463 


0.463 


-0.664 


-0.155 


0.108 


-0.00905 


0.0993 


2.8(-6) 




HF 


0.238 


0.238 


-0.325 


-0.0872 


0.064 


-0.00148 


0.0626 


0.00014 




DIRAC 


6.04 


6.04 


-6.71 


-7.15 


-1.77 


-1.21 


-2.98 


-0.00285 


( ( lO H 1 o 


PBEO 


4.58 


4.58 


-5.11 


-3.98 


0.0664 


-0.702 


-0.636 


-0.0412 


B3LYP 


4.86 


4.86 


-5.37 


-2.82 


1.54 


-0.708 


0.834 


0.047 




HF 


4.4 


A A 


— 4.yD 


— Z. 14 


1 7 
i. / 


— U.o ( 4 


I.OZ 


U.UUUoUl 




DIRAC 


16.3 


19.3 


-28.7 


-10.7 


-3.68 


-6.28 


-9.96 


-0.00402 


(HF)2 


PBEO 


13.3 


15.9 


-24.7 


-6.91 


-2.31 


-4.98 


-7.28 


-0.00026 


B3LYP 


14.2 


17 


-25.9 


-7 


-1.69 


-5.08 


-6.77 


-0.00223 




HF 


11 


13.1 


-22.2 


-4.17 


-2.19 


-3.61 


-5.8 


-2.12(-5) 




DIRAC 


11.8 


11.8 


-18.9 


-8.78 


-4.04 


-2.89 


-6.93 


-0.00268 


(NH3)2 


PBEO 


9.4 


9.4 


-15.7 


-5.49 


-2.36 


-2.06 


-4.41 


-0.000912 


B3LYP 


9.97 


9.97 


-16.4 


-4.9 


-1.31 


-2.11 


-3.43 


-0.000516 




HF 


8.69 


8.69 


-14.9 


-3.3 


-0.813 


-1.42 


-2.23 


-2.16(-6) 



^ -5e-06 

.g 

^§ -6e-06 




0.1 0.01 0.001 0.0001 le-05 le-06 

grid threshold/mH 

FIG. 1: Dependence of SEint on grid threshold for the Ne2 
dimer at ii = 6 ao for PBEO functional in aug-cc-pVQZ basis 
set. 



for all values of the grid threshold reported in Fig. 1. 
The calculations employed aug-cc-pVQZ basis set for 
neon dimer and aug-cc-pVTZ for the other systems. 
Dimer-centred basis set (DCBS) has been used through- 
out the calculations and the supermolecular interaction 
energies were counterpoise (CP)-corrected for the basis 
set superposition error (BSSE). 



B. Discussion 

From the results of Table I it is clear that the bifunc- 
tional approach converges to the same values of the inter- 
action energies as the conventional KS procedure, within 



excellent accuracy of below 0.003 % for the H-bondcd 
dimers, and 0.05 % for the van der Waals dimers. 

In general, as long as our procedure is convergent, it 
must converge to the same result as the standard KS ap- 
proach. This is because no extra constraints are imposed 
on the functional, and the finally optimized total density 
must be the same in both cases. However, the details of 
the convergence depend on several factors: 

• the particular functional, 

• the system under consideration, 

• the initial orthogonalization of monomer orbitals, 
and 

• the orthogonality forcing technique in the iteration 
process. 

The functional and system convergence dependence is ev- 
ident in Table I. The dependence on mutual orthognal- 
ization scheme and the manner it is forced in the iter- 
ation process have not been studied so far as only one 
approach has been adopted — they will be studied in 
the future in the context of actual applications. For in- 
stance, convergence problems may appear for complexes 
which undergo major redistribution of electron densities 
between monomers, such as in donor-acceptor interac- 
tions and also for those which are poorly described by a 
single determinant due to static correlation effects. 

The total interaction energy is composed of the HL en- 
ergy and the PB deformation energy. As pointed out in 
the previous Section, both and E^^ are uniquely 
defined, and are independent of the orthogonalization 
procedure. However, this is not the case for the first 
four terms of Eq. (33) in Table I. They describe the 
monomer deformation effects due to orthogonalization 
(cols. 3 and 4), and electrostatic and exchange-correlation 
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effects, arising between the monomers described with or- 
thogonalised occupied orbitals (cols. 5 and 6). Therefore 
they are not uniquely defined, and are strongly depen- 
dent on the orthogonalization scheme — they are not 
useful to interpret the interaction. 

The total interaction energies listed in the Table I are 
qualitatively correct only for polar molecules; for atoms 
and non-polar species they are erroneous due to the well 
known fact that the functionals: DIRAC, PBEO, B3LYP 
do not reproduce the dispersion component. The 
and components may be compared to similar terms 
at the HF level of theory (the last entry for every dimer) . 
Assuming that the functionals reproduce only the local 
correlation terms, but fail to recover dispersion contri- 
butions, the DFT E^ results should differ from the HF 
ones by a relatively small intranionomer correlation ef- 
fect. This is apparently not the case for DIRAC, for 
van der Waals complexes: Ne2 and (C2H4)2. Also PBEO 
shows attraction, albeit small, for Ne2, and seems to be 
not repulsive enough for (€2114)2. Such a behavior in- 
dicates that some residual dispersion terms are present. 
As to the B3LYP functional, it seems to be the most 
dispersion-free, since its values of E^ are the closest to 
the HF values. These results are in agreement with the 
observations that the amount of dispersion in a functional 
correlates to the steepness of the exchange enhancement 
factor: it becomes steeper when moving from DIRAC to 
B3LYP, and, consequently, the amount of dispersion in- 
cluded in these functionals is reduced (see Refs.^^"'^^) Fi- 
nally, for hydrogen-bonded complexes, all methods give 
qualitatively correct E^^, only PBEO and B3LYP are 
closer to the HF result than DIRAC. 

The above discussion suggests a useful application of 
the formalism presented in this paper. For approximate 
functionals the partitioning of the interaction energy into 
the E^ and components may serve as a diagnostic 
of the functional's adequacy in the intermolecular inter- 
action energy problems. Indeed, one can determine, what 
components and how efhciently are recovered by a tested 
fimctional. 



IV. SUMMARY AND OUTLOOK 

In this paper we provided a rigorous derivation of 
the supermolecular DFT interaction energy in terms of 



mutual monomer polarization via the Pauli blockade 
method. Numerical calculations for four model systems 
and three example functionals of different types have 
proved that the formalism leads to interaction energies 
rapidly converging to the supermolecular interaction en- 
ergies. The accuracy achieved is better than 0.1 % and 
appears to be limited only by the size of the grid. The 
accuracy is qualitatively similar for all three DFT func- 
tionals under investigation. Our formalism appears thus 
to be a viable and useful alternative of solving the KS 
equations in terms of separated-monomer orbitals rather 
than supermolecular orbitals. This fact has several im- 
portant practical implications. 

On the one hand, the presented formalism offers 
possibilities of using different functionals to describe 
the monomers and to describe the interaction. We 
have recently exploited this feature to design a novel 
DFT+dispersion approach.''^ On the other, it would be 
of interest to combine our approach for the DFT that is 
capable to reproduce the dispersion energy, e.g. of range- 
separated family of functionals.^ 

The bifunctional formulation provides a platform for 
deriving a choice of DFT treatments which use different 
functionals (or even theories) for different subsystems. 
The results may thus be of interest for those who use 
subsystem formulation in the context of embedding po- 
tentials. ^"■^^"'''■^ 

Our formalism may also be useful when interpret- 
ing different interaction energy decomposition schemes 
that have recently been proposed within the DFT ap- 
proach. ^^^^'''^'^ 
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